# Load Necessary Packages
library(glmnet)
library(foreign)
library(estimatr)
library(dplyr)
library(tidyr)
library(bcf)
library(dbarts)
library(texreg)
library(Hmisc)
library(grid)
library(gridExtra)
library(interplot)
library(margins)
library(car)
library(TOSTER)
library(interflex)
library(effectsize)

# Load Data
setwd("/users/josephphillips/Dropbox/COVID-19/Data/")
uk <- read.dta("GB_pooled_allwaves_v12.dta")

# Load GLMNet Function
run_LASSO <- function(data, yvar, pre){
  
  set.seed(2334)
  
  var_quo <- enquo(yvar)
  var_quo1 <- enquo(pre)
  
  ## select x variables 
  dat <- data %>%
    dplyr::select(university, agegroup, male, married, frequentchurch, profile_GOR, left, right, ideology1, highincidence, knowledge, nonwhite, polinterest,trust_healthgov,total_media_trust)
  
  pre <- data %>%
    dplyr::select(!!var_quo1) %>% pull() # Save as a vector instead of a data frame ot use model matrix
  
  ## select y-variable
  out1 <- data %>% 
    dplyr::select(!!var_quo) %>% pull() # save as a vector instead of a data-frame to use model.matrix()
  
  merged1 <- cbind(out1, dat, pre) %>% na.omit() ## merge and delete NA
  x1 <- model.matrix(out1 ~ ., data = merged1) ## make model matrix
  y1 <- merged1$out1 ## select y-variable from NA-deleted dataset
  mod1 <- cv.glmnet(x1, y1, alpha = 1, family = "gaussian") ## run lasso model
  coef <- coef(mod1) ## get coef  
  
  list(merged1, x1, y1, mod1, coef)
  
}

# Clean Variables
uk$treatment <- factor(ifelse(uk$splitvariable>=1 & uk$splitvariable<=3,"Contention",ifelse(uk$splitvariable>=4 & uk$splitvariable<=7,"Injunctive","Control")),levels=c("Control","Injunctive","Contention"))
uk$injunctive_treatment <- ifelse(uk$treatment=="Injunctive",1,0)
uk$contention_treatment <- ifelse(uk$treatment=="Contention",1,0)
uk$approve_Boris <- as.numeric(uk$approve_Boris)
uk$covid19_approve1 <- as.numeric(uk$covid19_approve1)
uk$polinterest <- as.numeric(uk$polinterest)
uk$covid19_exp_1 <- as.numeric(uk$covid19_exp_1)-1
uk$covid19_exp_3 <- as.numeric(uk$covid19_exp_3)-1
uk$covid19_exp_4 <- as.numeric(uk$covid19_exp_4)-1
uk$covid19_exp_5 <- as.numeric(uk$covid19_exp_5)-1
uk$covid19_exp_6 <- as.numeric(uk$covid19_exp_6)-1
uk$covid19_exp_7 <- as.numeric(uk$covid19_exp_7)-1
uk$agegroup <- ifelse(uk$agegroup==1,"age18-34",ifelse(uk$agegroup==2,"age35-44",ifelse(uk$agegroup==3,"age45-54",ifelse(uk$agegroup==4,"age55-64","age65+"))))
uk$total_media_trust <- uk$total_media_trust/2
uk$approve_Boris <- (uk$approve_Boris-1)/3
uk$therm_boris <- uk$therm_boris/100
uk$trust_healthgov <- (uk$trust_healthgov-1)/3
uk$conspiracy_general <- 1-(uk$conspiracy_general-1)/4
uk$total_media_trust <- (uk$total_media_trust-1)/3
uk$knowledge <- uk$knowledge/5
uk$lagdv <- as.numeric(uk$covid19_vaccine)/5
uk$dv_missing <- ifelse(uk$W3_covid19_vaccine=="N/A, already received the vaccine",NA,7-as.numeric(uk$W3_covid19_vaccine))
uk$dv_verylikely <- ifelse(uk$W3_covid19_vaccine=="N/A, already received the vaccine",6,7-as.numeric(uk$W3_covid19_vaccine))
uk$dv_mostlikely <- ifelse(uk$W3_covid19_vaccine=="N/A, already received the vaccine",7,7-as.numeric(uk$W3_covid19_vaccine))
uk$age65 <- ifelse(uk$agegroup=="age65+",1,0)
uk$vaccine_norms_guess <- ifelse(uk$vaccine_norms_guess>100,NA,uk$vaccine_norms_guess/100)
uk$underestimated_norms <- ifelse(uk$vaccine_norms_guess<0.72,1,0)
uk$overestimated_norms <- ifelse(uk$vaccine_norms_guess>0.90,1,0)
uk$factcheck_treatment <- 2-uk$factcheck_treatment
uk$W3factcheck_treatment <- 2-uk$W3factcheck_treatment
uk$contention <- (5-as.numeric(uk$W2_vaccine_contention))/4
uk$contention_accuracy <- ifelse(uk$W2_vaccine_contention=="Strongly agree" | uk$W2_vaccine_contention=="Somewhat agree",1,0)
uk$agegroup <- ifelse(uk$age>=18 & uk$age<=34,"18 to 34",ifelse(uk$age>=35 & uk$age<=49,"35 to 49",ifelse(uk$age>=50 & uk$age<=64,"50 to 64","65 and older")))

# Lasso
df <- uk
lasso_missing <- run_LASSO(df,dv_missing,lagdv) ## Health Trust, Pre
lasso_verylikely <- run_LASSO(df,dv_verylikely,lagdv) ## Health Trust, Pre
lasso_mostlikely <- run_LASSO(df,dv_mostlikely,lagdv) ## Age 65+, Health Trust, Pre

# Table S2
table(uk$treatment,uk$university)
chisq.test(uk$treatment,uk$university)
table(uk$treatment,uk$agegroup)
uk$age1834 <- ifelse(uk$agegroup=="18 to 34",1,0)
uk$age3549 <- ifelse(uk$agegroup=="35 to 49",1,0)
uk$age5064 <- ifelse(uk$agegroup=="50 to 64",1,0)
uk$age65 <- ifelse(uk$agegroup=="65 and older",1,0)
chisq.test(uk$treatment,uk$age1834)
chisq.test(uk$treatment,uk$age3549)
chisq.test(uk$treatment,uk$age5064)
chisq.test(uk$treatment,uk$age65)
table(uk$treatment,uk$male)
chisq.test(uk$treatment,uk$male)
table(uk$treatment,uk$married)
chisq.test(uk$treatment,uk$married)
table(uk$treatment,uk$frequentchurch)
chisq.test(uk$treatment,uk$frequentchurch)
table(uk$treatment,uk$profile_GOR)
uk$northeast <- ifelse(uk$profile_GOR=="North East",1,0)
uk$northwest <- ifelse(uk$profile_GOR=="North West",1,0)
uk$yorkshire <- ifelse(uk$profile_GOR=="Yorkshire and the Humber",1,0)
uk$eastmidlands <- ifelse(uk$profile_GOR=="East Midlands",1,0)
uk$westmidlands <- ifelse(uk$profile_GOR=="West Midlands",1,0)
uk$east <- ifelse(uk$profile_GOR=="East of England",1,0)
uk$london <- ifelse(uk$profile_GOR=="London",1,0)
uk$southeast <- ifelse(uk$profile_GOR=="South East",1,0)
uk$southwest <- ifelse(uk$profile_GOR=="South West",1,0)
uk$wales <- ifelse(uk$profile_GOR=="Wales",1,0)
uk$scotland <- ifelse(uk$profile_GOR=="Scotland",1,0)
chisq.test(uk$treatment,uk$northeast)
chisq.test(uk$treatment,uk$northwest)
chisq.test(uk$treatment,uk$yorkshire)
chisq.test(uk$treatment,uk$eastmidlands)
chisq.test(uk$treatment,uk$westmidlands)
chisq.test(uk$treatment,uk$east)
chisq.test(uk$treatment,uk$london)
chisq.test(uk$treatment,uk$southeast)
chisq.test(uk$treatment,uk$southwest)
chisq.test(uk$treatment,uk$wales)
chisq.test(uk$treatment,uk$scotland)
table(uk$treatment,uk$right)
chisq.test(uk$treatment,uk$right)
table(uk$treatment,uk$left)
chisq.test(uk$treatment,uk$left)
mean(subset(uk,treatment=="Control")$ideology1,na.rm=T)
mean(subset(uk,treatment=="Injunctive")$ideology1,na.rm=T)
mean(subset(uk,treatment=="Contention")$ideology1,na.rm=T)
summary(aov(ideology1~treatment,data=uk))
table(uk$treatment,uk$highincidence)
chisq.test(uk$treatment,uk$highincidence)
mean(subset(uk,treatment=="Control")$knowledge*5,na.rm=T)
mean(subset(uk,treatment=="Injunctive")$knowledge*5,na.rm=T)
mean(subset(uk,treatment=="Contention")$knowledge*5,na.rm=T)
summary(aov(knowledge~treatment,data=uk))
table(uk$treatment,uk$nonwhite)
chisq.test(uk$treatment,uk$nonwhite)
mean(subset(uk,treatment=="Control")$polinterest,na.rm=T)
mean(subset(uk,treatment=="Injunctive")$polinterest,na.rm=T)
mean(subset(uk,treatment=="Contention")$polinterest,na.rm=T)
summary(aov(polinterest~treatment,data=uk))
mean(subset(uk,treatment=="Control")$trust_healthgov*3,na.rm=T)
mean(subset(uk,treatment=="Injunctive")$trust_healthgov*3,na.rm=T)
mean(subset(uk,treatment=="Contention")$trust_healthgov*3,na.rm=T)
summary(aov(trust_healthgov~treatment,data=uk))
mean(subset(uk,treatment=="Control")$total_media_trust*3,na.rm=T)
mean(subset(uk,treatment=="Injunctive")$total_media_trust*3,na.rm=T)
mean(subset(uk,treatment=="Contention")$total_media_trust*3,na.rm=T)
summary(aov(total_media_trust~treatment,data=uk))
mean(subset(uk,treatment=="Control")$lagdv,na.rm=T)
mean(subset(uk,treatment=="Injunctive")$lagdv,na.rm=T)
mean(subset(uk,treatment=="Contention")$lagdv,na.rm=T)
summary(aov(lagdv~treatment,data=uk))

# Table S4/Figure 1
Reg.h1a.verylikely.uk.robust <- lm_robust(dv_verylikely~treatment+trust_healthgov+lagdv,data=uk)
uk1 <- effectsize(Reg.h1a.verylikely.uk.robust)
Reg.h1a.mostlikely.uk.robust <- lm_robust(dv_mostlikely~treatment+trust_healthgov+age65+lagdv,data=uk)
Reg.h1a.missing.uk.robust <- lm_robust(dv_missing~treatment+trust_healthgov+lagdv,data=uk)

# Table S7/Figure 2
Reg.h1b.verylikely.uk.robust <- lm_robust(dv_verylikely~(injunctive_treatment*underestimated_norms)+(injunctive_treatment*overestimated_norms)+contention_treatment+trust_healthgov+lagdv,data=uk)
Reg.h1b.uk.under <- summary(margins(Reg.h1b.verylikely.uk.robust,at=list(underestimated_norms=1)))
Reg.h1b.uk.over <- summary(margins(Reg.h1b.verylikely.uk.robust,at=list(overestimated_norms=1)))
Reg.h1b.mostlikely.uk.robust <- lm_robust(dv_mostlikely~(injunctive_treatment*underestimated_norms)+(injunctive_treatment*overestimated_norms)+contention_treatment+trust_healthgov+age65+lagdv,data=uk)
Reg.h1b.missing.uk.robust <- lm_robust(dv_missing~(injunctive_treatment*underestimated_norms)+(injunctive_treatment*overestimated_norms)+contention_treatment+trust_healthgov+lagdv,data=uk)

# Table S9/Figure 3
Reg.rq1.verylikely.uk.robust <- lm_robust(dv_verylikely~+(contention_treatment*contention_accuracy)+injunctive_treatment+trust_healthgov+lagdv,data=uk)
Reg.rq1.uk.accurate <- summary(margins(Reg.rq1.verylikely.uk.robust,at=list(contention_accuracy=1)))
Reg.rq1.mostlikely.uk.robust <- lm_robust(dv_mostlikely~+(contention_treatment*contention_accuracy)+injunctive_treatment+trust_healthgov+age65+lagdv,data=uk)
Reg.rq1.missing.uk.robust <- lm_robust(dv_missing~+(contention_treatment*contention_accuracy)+injunctive_treatment+trust_healthgov+lagdv,data=uk)

# Table S10
Reg.rq2.verylikely.uk.robust <- lm_robust(dv_verylikely~(injunctive_treatment*factcheck_treatment)+(injunctive_treatment*W3factcheck_treatment)+(contention_treatment*factcheck_treatment)+(contention_treatment*W3factcheck_treatment)+(contention_treatment*factcheck_treatment)+(contention_treatment*W3factcheck_treatment)+trust_healthgov+lagdv,data=uk)
linearHypothesis(Reg.rq2.verylikely.uk.robust,c("factcheck_treatment=0","W3factcheck_treatment=0","injunctive_treatment:factcheck_treatment=0","injunctive_treatment:W3factcheck_treatment=0","factcheck_treatment:contention_treatment=0","W3factcheck_treatment:contention_treatment=0"),test="F") # p=.996, fail to reject null
Reg.rq2.mostlikely.uk.robust <- lm_robust(dv_mostlikely~(injunctive_treatment*factcheck_treatment)+(injunctive_treatment*W3factcheck_treatment)+(contention_treatment*factcheck_treatment)+(contention_treatment*W3factcheck_treatment)+trust_healthgov+age65+lagdv,data=uk)
linearHypothesis(Reg.rq2.mostlikely.uk.robust,c("factcheck_treatment=0","W3factcheck_treatment=0","injunctive_treatment:factcheck_treatment=0","injunctive_treatment:W3factcheck_treatment=0","factcheck_treatment:contention_treatment=0","W3factcheck_treatment:contention_treatment=0"),test="F") # p=.992, fail to reject null
Reg.rq2.missing.uk.robust <- lm_robust(dv_missing~(injunctive_treatment*factcheck_treatment)+(injunctive_treatment*W3factcheck_treatment)+(contention_treatment*factcheck_treatment)+(contention_treatment*W3factcheck_treatment)+trust_healthgov+lagdv,data=uk)
linearHypothesis(Reg.rq2.missing.uk.robust,c("factcheck_treatment=0","W3factcheck_treatment=0","injunctive_treatment:factcheck_treatment=0","injunctive_treatment:W3factcheck_treatment=0","factcheck_treatment:contention_treatment=0","W3factcheck_treatment:contention_treatment=0"),test="F") # p=.994, fail to reject null

# Table S12
Reg.h1a.verylikely.nocov.uk <- lm_robust(dv_verylikely~treatment,data=uk)
Reg.h1a.mostlikely.nocov.uk <- lm_robust(dv_mostlikely~treatment,data=uk)
Reg.h1a.missing.nocov.uk <- lm_robust(dv_missing~treatment,data=uk)

# Table S15
Reg.h1b.verylikely.uk.robust.nocov <- lm_robust(dv_verylikely~(injunctive_treatment*underestimated_norms)+(injunctive_treatment*overestimated_norms),data=uk)
Reg.h1b.mostlikely.uk.robust.nocov <- lm_robust(dv_verylikely~(injunctive_treatment*underestimated_norms)+(injunctive_treatment*overestimated_norms),data=uk)
Reg.h1b.missing.uk.robust.nocov <- lm_robust(dv_verylikely~(injunctive_treatment*underestimated_norms)+(injunctive_treatment*overestimated_norms),data=uk)

# Table S17
Reg.rq1.verylikely.uk.robust.nocov <- lm_robust(dv_verylikely~+(contention_treatment*contention_accuracy),data=uk)
Reg.rq1.mostlikely.uk.robust.nocov <- lm_robust(dv_mostlikely~+(contention_treatment*contention_accuracy),data=uk)
Reg.rq1.missing.uk.robust.nocov <- lm_robust(dv_missing~+(contention_treatment*contention_accuracy),data=uk)

# Table S18
Reg.rq2.verylikely.nocov.uk.robust <- lm_robust(dv_verylikely~(injunctive_treatment*factcheck_treatment)+(injunctive_treatment*W3factcheck_treatment)+(contention_treatment*factcheck_treatment)+(contention_treatment*W3factcheck_treatment)+(contention_treatment*factcheck_treatment)+(contention_treatment*W3factcheck_treatment),data=uk)
linearHypothesis(Reg.rq2.verylikely.nocov.uk.robust,c("factcheck_treatment=0","W3factcheck_treatment=0","injunctive_treatment:factcheck_treatment=0","injunctive_treatment:W3factcheck_treatment=0","factcheck_treatment:contention_treatment=0","W3factcheck_treatment:contention_treatment=0"),test="F") # p=.954, fail to reject null
Reg.rq2.mostlikely.nocov.uk.robust <- lm_robust(dv_mostlikely~(injunctive_treatment*factcheck_treatment)+(injunctive_treatment*W3factcheck_treatment)+(contention_treatment*factcheck_treatment)+(contention_treatment*W3factcheck_treatment),data=uk)
linearHypothesis(Reg.rq2.mostlikely.nocov.uk.robust,c("factcheck_treatment=0","W3factcheck_treatment=0","injunctive_treatment:factcheck_treatment=0","injunctive_treatment:W3factcheck_treatment=0","factcheck_treatment:contention_treatment=0","W3factcheck_treatment:contention_treatment=0"),test="F") # p=.997, fail to reject null
Reg.rq2.missing.nocov.uk.robust <- lm_robust(dv_missing~(injunctive_treatment*factcheck_treatment)+(injunctive_treatment*W3factcheck_treatment)+(contention_treatment*factcheck_treatment)+(contention_treatment*W3factcheck_treatment),data=uk)
linearHypothesis(Reg.rq2.missing.nocov.uk.robust,c("factcheck_treatment=0","W3factcheck_treatment=0","injunctive_treatment:factcheck_treatment=0","injunctive_treatment:W3factcheck_treatment=0","factcheck_treatment:contention_treatment=0","W3factcheck_treatment:contention_treatment=0"),test="F") # p=.704, fail to reject null

# Figure S1, GB
Reg.h1b.verylikely.uk.flex <- interflex("kernel",data=uk,Y="dv_verylikely",D="injunctive_treatment",X="vaccine_norms_guess",Z=c("contention_treatment","trust_healthgov","lagdv"),na.rm=T,CI=T)
